function [A, B, C, D, E] = CMTF( X, Y, A, B, C, D, E, mu, sigma, Err, maxIter)
    
    % default
    % X - 4 dimension & Y - 3 dimension
    
    % X (A, B, C, D)
    % input: A, B, C, D, X non-empty matrix or tensor
    
    % Y (A, B, E)
    % input: A, B, E, Y non-empty matrix or tensor
    
    dimType = 1;
    
    % X - 3 dimension & Y - 3 dimension
    
    % X (A, B, C)
    % input: A, B, C, X non-empty matrix or tensor
    %        D = [] 
    
    % Y (A, B, E)
    % input: A, B, E, Y non-empty matrix or tensor
    
    if(isempty(D))
        dimeType = 2;
    end
    
    tempSize = size(A);
    R = tempSize(2);
    I = tempSize(1);
    tempSize = size(B);
    J = tempSize(1);
    tempSize = size(C);
    K = tempSize(1);
    if(dimeType == 1)
        tempSize = size(D);
        M = tempSize(1);
    else
        M = 0;
    end
    tempSize = size(E);
    L = tempSize(1);
    
    
    
    % change to err for convergence
    iter = 0;
    err = Err;
    e1 = 0;
    e2 = 0;
    
    while(iter < maxIter)
        if(dimType == 1)
            
            % calculate e1
            for i = 1:I
                for j = 1:J
                    for k = 1:K
                        for m = 1:M
                            elewiseProduct = sum((A(i,:).*B(j,:)).*(C(k,:).*D(m,:)));
                            e1 = e1 + (X(i,j,k,m) - elewiseProduct) * (X(i,j,k,m) - elewiseProduct);
                        end
                    end
                end
            end
            e1 = (1/sigma) * e1;
            % calculate e2
            for i = 1:I
                for j = 1:J
                    for l = 1:L
                        elewiseProduct = sum((A(i,:).*B(j,:)).*E(l,:));
                        e2 = e2 + (Y(i,j,l) - elewiseProduct) * (Y(i,j,l) - elewiseProduct);
                    end
                end
            end
            e2 = (1/sigma) * e2;
            
            % update A
            for i = 1 : I
                for r = 1 : R
                    part1 = sum(B(:,r))*sum(C(:,r))*sum(D(:,r));
                    part2 = sum(B(:,r))*sum(E(:,r));
                    A(i,r) = A(i,r) + mu*(e1*part1 + e2*part2);
                end
            end
            % update B
            for j = 1 : J
                for r = 1 : R
                    part1 = sum(A(:,r))*sum(C(:,r))*sum(D(:,r));
                    part2 = sum(A(:,r))*sum(E(:,r));
                    B(j,r) = B(j,r) + mu*(e1*part1 + e2*part2);
                end
            end
            % update C
            for k = 1 : K
                for r = 1 : R
                    part1 = sum(A(:,r))*sum(B(:,r));
                    C(k,r) = C(k,r) + mu*(e1*part1);
                end
            end
            % update D
            for m = 1 : M
                for r = 1 : R
                    part1 = sum(A(:,r))*sum(B(:,r));
                    D(m,r) = D(m,r) + mu*(e1*part1);
                end
            end
            % update E
            for l = 1 : L
                for r = 1 : R
                    part1 = sum(A(:,r))*sum(B(:,r));
                    E(l,r) = E(l,r) + mu*(e2*part1);
                end
            end
        else
            
            % calculate e1
            for i = 1:I
                for j = 1:J
                    for k = 1:K
                        elewiseProduct = sum((A(i,:).*B(j,:)).*C(k,:));
                        e1 = e1 + (X(i,j,k) - elewiseProduct) * (X(i,j,k) - elewiseProduct);
                    end
                end
            end
            e1 = (1/sigma) * e1;
            % calculate e2
            for i = 1:I
                for j = 1:J
                    for l = 1:L
                        elewiseProduct = sum((A(i,:).*B(j,:)).*E(l,:));
                        e2 = e2 + (Y(i,j,l) - elewiseProduct) * (Y(i,j,l) - elewiseProduct);
                    end
                end
            end
            e2 = (1/sigma) * e2;
            
            % update A
            for i = 1 : I
                for r = 1 : R
                    part1 = sum(B(:,r))*sum(C(:,r));
                    part2 = sum(B(:,r))*sum(E(:,r));
                    A(i,r) = A(i,r) + mu*(e1*part1 + e2*part2);
                end
            end
            % update B
            for j = 1 : J
                for r = 1 : R
                    part1 = sum(A(:,r))*sum(C(:,r));
                    part2 = sum(A(:,r))*sum(E(:,r));
                    B(j,r) = B(j,r) + mu*(e1*part1 + e2*part2);
                end
            end
            % update C
            for k = 1 : K
                for r = 1 : R
                    part1 = sum(A(:,r))*sum(B(:,r));
                    C(k,r) = C(k,r) + mu*(e1*part1);
                end
            end
            % update E
            for l = 1 : L
                for r = 1 : R
                    part1 = sum(A(:,r))*sum(B(:,r));
                    E(l,r) = E(l,r) + mu*(e2*part1);
                end
            end
        end
    end
end